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Preface 

This  thesis  is  a  continuation  of  previous  work  done 
at  the  Air  Force  Institute  of  Technology  on  Monte  Carlo  tech¬ 
niques  of  reliability  prediction.  It  is  hoped  that  the 
method  proposed  here  may  find  use  in  reliability  analysis  of 
complex  systems  and/or  in  situations  where  limited  failure 
test  data  is  available  to  the  reliability  engineer. 

I  wish  to  thank  my  thesis  advisor.  Professor  Albert  H. 
Moore,  for  suggesting  the  topic  and  giving  guidance  and 
encouragement  throughout  this  project. 

I  also  wish  to  express  my  gratitude  to  Dr.  H.  Leon 
Harter  who  gave  valuable  advise  during  this  project.  The 
FORTRAN  programs  in  Appendix  A  and  Appendix  C  are  revisions 
of  programs  which  he  originally  wrote  and  kindly  supplied 
to  me , 
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Abstract 

This  paper  develops  a  Monte  Carlo  technique  which,  with 
a  digital  computer,  determines  confidence  limits  for  system 
reliability  of  complex  systems  containing  components  charac¬ 
terized  by  the  Weibull  distribution.  The  component  distribu¬ 
tion  shape  and  scale  parameters  are  estimated  by  the  method 
of  maximum  likelihood  from  component  failure  times  while  the 
location  parameter  is  assumed  known.  The  asymptotic  distribu¬ 
tion  of  these  maximum  likelihood  estimators  and  a  Monte  Carlo 
simulation  are  used  to  determine  confidence  limits  on  system 
reliability.  As  an  example,  confidence  limits  are  calculated 
for  two  systems  of  up  to  eight  components  in  combinations  of 
series  and  parallel  configurations  using  99,  499,  999,  and 
2999  simulations.  Accuracy  of  the  confidence  limits  is  found 
to  be  satisfactory  after  being  checked  by  a  method  using  a 
double  Monte  Carlo  technique  which  assumes  values  for 
parameters  of  the  component  distributions  and  generates 
component  failure  times  to  be  used  in  the  proposed  Monte 
Carlo  technique.  Central  processing  time  for  this  technique 
on  the  CDC  6600  digital  computer  is  found  to  be  less  than  two 
minutes  for  the  examples  given.  There  is  no  limit  to  the 
number  of  components  or  the  configurations  of  the  components 
for  the  systems  to  which  the  method  can  be  applied.  The 
FORTRAN  IV  computer  programs  used  are  given  along  with  flow¬ 
charts  to  facilitate  reading  the  programs. 
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I .  Introduction 

The  Problem 

The  purpose  of  this  thesis  is  to  determine  whether 
accurate  interval  estimation  of  reliability  for  complicated 
systems,  containing  dissimilar  components  characterized  by 
the  Weibull  distribution,  can  be  obtained  by  using  maximum 
likelihood  estimates  and  a  Monte  Carlo  simulation. 

Significance .  Systems  containing  many  components  become 
increasingly  complex  and  costly.  Complicated  systems  costing 
thousands  of  dollars  can  become  inoperative  due  to  the 
failure  of  a  single  connection.  Because  of  this,  it  is 
important  to  determine  system  reliability  as  a  function  of 
the  component  reliabilities  which  make  up  the  system.  This 
can  be  done  fairly  easily  for  some  few  types  of  combinations 
of  components  which  are  from  certain  parent  populations  such 
as  an  electrical  circuit  made  up  of  components  characterized 
by  the  exponential  distribution.  Khen  the  components  are 
dissimilar,  no  simple  analytical  expression  can  be  used  to 
compute  system  reliability. 

In  addition,  predictions  of  reliability  should  be  based 
on  test  data.  However,  as  system  reliability  requirements 
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increase,  testing  a  system  to  destruction  can  become  expen¬ 
sive,  both  in  time  needed  and  money  spent. 

A  method  that  circumvents  this  problem  uses  limited 
test  data  to  calculate  maximum  likelihood  estimates  of  the 
parameters  of  the  failure  distribution.  These  estimates  are 
used,  along  with  a  random  number  generator,  in  a  Monte  Carlo 
technique  which  uses  the  asymptotic  distribution  of  the 
maximum  likelihood  estimators  to  generate  more  point  samples 
of  system  reliability.  These  point  samples  are  used  to 
determine  an  interval  estimate  of  system  reliability  with 
the  associated  confidence  level. 

Definitions .  Reliability  is  the  probability  that  a 


system  is  still  operating  at  time  t.  If  T  is  the  time  to 
failure  or  life  length  of  a  system  or  component,  the  relia¬ 
bility  at  time  t  or  R(t)  is  given  by  R(t)  =  P (T  >  t) ,  where 
P  means  "probability  of". 

The  Weibull  density  function  is  defined  as 


f (t;C,e,K) 


K(t-C)K-1 

0K 


(11 


0,  K  >  0,  C  <  t 


=  0  elsewhere 


where  0  =  the  scale  parameter,  K  =  the  shape  parameter,  and 
C  *=  location  parameter.  The  scale  parameter  affects  the 
dispersion  of  the  random  variable  t  about  its  mean.  The 
shape  parameter  determines  whether  the  hazard  function  is 
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R(t)  =  exp  J—-]  |  (2) 

Pseudo-random  numbers  are  generated  from  a  series  of 
numbers  which  repeats  itself  with  long  periods  between 
repetitions,  "hese  numbers  are  suitably  random  if  the 
repetition  periods  are  greater  than  the  number  of  digits 
desired  and  if  there  is  a  small  enough  correlation  between 
digits  (Ref  15:257). 

The  maximum  likelihood  estimator  (MLE)  is  the  estimator 
which  maximizes  the  likelihood  function  with  respect  to  the 
parameter  being  estimated.  The  likelihood  function  is  the 
joint  density  of  a  sample  of  n  random  variables  and  is 
expressed  by 

n 

L(x1,x2.  .xn;9)  =  fCx^e)  (3) 

where  L(x.:8)  =  the  likelihood  function. 

i 

A 

If  it  is  desired  to  determine  0  (the  MLE  for  the 
parameter  0),  normally  it  is  only  necessary  to  take  the 
partial  derivative  of  the  likelihood  function  with  respect 
to  0,  set  it  equal  to  zero  and  solve  for  0  which  then  becomes 

A 

0.  By  using  the  MLE  for  each  parameter  and  the  expression 
for  reliability  for  a  given  failure  density,  it  is  possible 
to  arrive  at  a  point  estimate  for  reliability. 
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The  method  in  this  paper  generates  intervals  from  a 
random  sample  of  system  reliabilities  and  this  is  called  a 
random  interval.  If  reliability  falls  within  any  given 
interval,  there  is  associated  with  the  interval  a  confidence 
level.  For  example,  if  it  is  said  the  reliability  of  a 
system  lies  between  0.6  and  0.8  at  the  95%  confidence  level, 
it  is  meant  P(0.6  <  R(t)  <  0.8)  *  .95  or  the  random  interval 
0.6  -  0.8  would  contain  the  true  system  reliability  95%  of 
the  time  if  the  system  were  tested  many  times.  The  interval 
is  called  a  confidence  interval,  which  has  0.6  as  a  lower 
confidence  limit  and  0.8  as  an  upper  confidence  limit. 

A  Monte  Carlo  simulation  is  one  in  which  a  random  value 
or  number  is  picked  from  all  possible  numbers  that  can  be 
described  by  a  specific  density  function. 


Analysis 

Assumptions .  This  thesis  will  concern  itself  with 
systems  which  have  components  described  by  the  Weibull 
density  function.  It  is  assumed  that  the  components  have 
previously  been  determined  to  be  Weibull  or  that  the  Weibull 
density  function  adequately  models  the  components  in  the 
system.  This  paper  does  not  provide  a  method  to  test  for 
the  Weibull  density  function  though  these  tests  do  exist  and 
are  described  in  other  references.  One  of  the  easiest  ways 
to  test  for  the  Weibull  distribution  is  to  use  Weibull  graph 
paper  (Ref  5:170-173).  Once  it  is  determined  or  suspected 
the  components  are  Weibull,  the  methods  in  this  paper  can  be 
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used  to  estimate  the  parameters  of  the  Weibull  distribution 
a.d  then  to  generate  interval  estimates  of  system  reliability. 
The  Weibull  distribution  is  rather  versatile  in  that, 
depending  on  the  value  of  the  shape  parameters,  the  Weibull 
distribution  can  be  used  to  model  components  with  decreasing, 
increasing  or  constant  failure  rates  or  hazard  functions. 

A  special  case  of  the  Weibull  distribution  is  the  exponential 
distribution  which  occurs  when  the  shape  parameter,  K,  is 
equal  to  one. 

It  is  assumed  that  the  components  of  the  system  being 
analyzed  fail  independently.  This  allows  use  of  standard 
formulas  for  components  connected  in  series  or  in  parallel. 

For  instance,  if  a  system  is  composed  of  two  components 
connected  in  series,  the  ?„*eliability  can  be  expressed  by 

Rs(t)  =  lRj(t)][R2Ct)]  (4) 

where  R ^  ( t )  =  reliability  of  component  number  one  and, 

R^Ct)  =  reliability  ox  exponent  number  two. 

Likewise,  for  a  system  made  up  of  two  components  in 
parallel  the  reliability  can  be  expressed  by 

Rs(t)  »  i  -  u  -  ^ctmi  -  R2(t)j 

*  i  -  QjQ2  (S) 
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where  =  1  -  R .  (t)  . 

More  complex  systems  can  bo  reduced  to  combinations  of 
series  and/or  parallel  configurations  by  use  of  Bayes* 
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theorem  or  the  Boolean  Disjunctive  theorem.  In  other  cases 
where  the  system  reliability  dependence  on  component  relia¬ 
bility  is  known  and  can  be  expressed  by  other  equations,  the 
independent  failure  assumption  is  not  necessary. 

It  is  assumed  that  the  location  parameter  is  known  or 
can  be  set  equal  to  zero  for  all  Weibull  component  distribu¬ 
tions  in  this  paper. 

Standards .  The  only  true  criterion  that  is  applied  to 
test  the  acceptability  of  the  method  in  this  paper  is  analysis 
of  the  result.  Does  the  computed  confidence  interval  contain 
the  actual  value  of  system  reliability  the  required  number  of 
times?  To  check  the  validity  of  the  method,  the  values  of 
all  parameters  were  assumed  for  each  of  the  Weibull  distribu¬ 
tions  of  the  components.  Pseudo-random  component  failure 
times  were  generated  for  each  of  the  components  and  the 
method  in  this  paper  applied  to  generate  many  confidence 
intervals.  Since  the  true  system  reliability  was  known,  a 
check  was  made  to  see  if  the  true  reliability  was  in  the 
calculated  intervals  the  required  number  of  times.  This 
process  is  described  in  Chapter  II. 

For  the  estimation  of  the  parameters  from  the  generated 
component  failure  times,  variation  from  the  true  value  of 
the  parameter  in  the  second  digit  was  deemed  acceptable. 

Background .  Levy  and  Moore  developed  a  method  for 
determining  system  reliability  confidence  limits  using  a 
Monte  Carlo  method  where  the  component  failure  distri outions 
(exponential,  normal,  lognormal.  Gamma,  or  Weibull)  with 
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location  and  shape  parameters  were  known  (Ref  10).  Lutton 
extended  this  method  by  using  the  asymptotic  distribution  of 
the  maximum  likelihood  estimators  of,  among  others,  the 
Weibull  distribution,  location  parameter  known  (Ref  11). 

Harter  and  Moore  developed  two  procedures  relevant  to 
this  paperi  They  used  an  iterative  procedure  to  determine 
maximum  likelihood  estimates  of  the  parameters  of  the  Weibull 
distribution  from  complete  and  censored  samples  (Ref  8). 

They  also  developed  a  method  to  determine  the  asymptotic 
variances  and  covariances  of  the  maximum  likelihood  estima¬ 
tors  of  the  parameters  of  the  Weibull  distribution  from 
samples  that  were  censored  or  complete  (Ref  7).  This  paper 
presents  a  technique  in  which  the  sample  of  failure  times 
for  each  component  is  used,  along  with  the  asymptotic  distri¬ 
bution  of  the  Weibull  parameters,  to  generate  a  confidence 
interval  with  the  associated  confidence  level. 

Approach .  The  method  and  associated  computer  program 
in  this  paper: 

1.  Calculate  the  MLE  of  the  component  shape  and  scale 
parameters  from  the  component  failure  times- 

2.  Determine  the  variance-covariance  matrix  of  the 
asymptotic  normal  distribution  of  the  shape  and 
scale  parameters. 

3.  Generate  a  sample  of  the  component  shape  and  scale 
parameters . 

4.  Use  sample  shape  and  scale  parameters,  along  with 
known  location  parameter,  to  calculate  a  point 
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estimate  of  component  reliability  for  each 
component . 

5.  Calculate  a  point  sample  of  system  reliability  from 
the  point  samples  of  component  reliabilities. 

6.  Repeat  steps  3-5  to  obtain  many  samples  of  system 
reliability. 

7.  Order  the  system  reliability  samples  to  obtain  an 
interval  sample  of  system  reliability  from  which 
the  confidence  limits  are  determined  at  a  given 
confidence  level. 
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1 1 .  Procedure 

Estimation  of  Component  Parameters 

The  computer  program  for  finding  the  MLE  of  the  para¬ 
meters  of  the  Weibull  distribution  from  complete  or  censored 
samples  is  given  in  Appendix  A.  The  procedure,  as  developed 
by  Harter  and  Moore,  uses  the  natural  logarithm  of  the 
likelihood  function  of  the  remaining  order  statistics  which 
is 

«  in  n!  -  Jln(n-ra)!  -  fcn  t!  +  (m-r)  (Jin  K  -  K  Jin  0)  + 


where  (m-r)  is  the  number  of  order  statistics  x  , ,  x  x 

r+1  r+2  m 

and  n  is  the  size  of  the  sample  before  censoring. 

The  partial  derivatives  of  with  respect  to  each  of 
the  three  parameters  arc  set  equal  to  zero  and  the  three 
resulting  equations  solved  simultaneously.  The  iterative 
procedure  for  the  simultaneous  solution  of  the  equations  uses 
an  initial  estimate  of  the  parameter  (which  is  fed  into  the 
computer  along  with  the  failure  times)  and  the  rule  of  false 
position  to  determine  a  new  value  of  the  parameter  from  the 
appropriate  likelihood  equation  into  which  the  latest  esti¬ 
mates  (or  a  known  value)  of  the  other  twu  parameters  have  been 
substituted  (Ref  8:641). 
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so  that  the  density  function  of  2  is  made  up  of  independent 
variates  of  the  normal  distribution  with  mean  zero  and 
variance  1  (Ref  14).  Thus,  to  obtain  a  sample  of  the  joint 
distribution  of  the  estimators,  the  following  formula  is 
used : 


Y' 


V1/2Z  +  Y 


(ID 


where  Y1  =  2  X  1  vector  of  the  biased  random  samples  of  the 
estimates  of  the  shape  and  scale  parameters. 

V  *  2  X  2  variance-covariance  matrix. 

Z  *=  2  X  1  vector  of  pseudo-random  numbers  from  the 
standard  normal  distribution. 

Y  =  2  X  1  vector  of  the  maximum  likelihood  estimates 

of  the  two  parameters  (as  found  in  the  previous 
section)  . 


or 

(12) 

In  formula  (12)  the  MLF,  of  9  and  K  are  used  as  the  marginal 
means  of  the  bivariate  normal  distribution  of  the  parameter 
estimates. 

The  random  samples  of  the  component  parameter  estimates 
were  substituted  into  the  expression  for  reliability  for 
each  component  to  get  random  values  of  component  reliabilities. 
Tht  were  then  substituted  into  the  expression  for  system 
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reliability,  one  value  for  each  component,  to  generate 

random  values  of  system  reliability. 

The  program  which  generated  the  system  reliabilities  is 

given  in  Appendix  B,  along  with  a  flow  chart  of  the  program. 

1/2 

A  detailed  explanation  of  the  evaluation  of  V  is  also 
given  in  Appendix  B. 


The  Variance -Covariance  Matrix 

The  natural  logarithm  of  the  likelihood  function  for  a 

sample  of  size  n  from  a  Weibull  population  where  the  lowest 

r  and  the  highest  n-m  sample  values  have  been  censored  is 

given  in  Ref  7.  The  elements  of  the  information  matrix  are 

found  by  taking  the  limits,  as  n  00 ,  of  the  negatives  of 

the  expected  values  of  the  second  partial  derivatives  of  the 

likelihood  function  with  respect  to  each  of  the  parameters. 

Let  k.  .  be  the  elements  of  a  matrix  B,  formed  from  the 

elements  of  the  information  matrix  (s^,  *  s21*  ant*  S22^’ 

where  E[...]  is  the  conditional  expectation  given  zy+j  anc* 

z  .  Then 
m 


— —  «  lim  pr  E  f  -  f  — 
n  r  In 


Kp  ♦  K(K+l)[r(2;z*)  -  T(2;^+1)]  ♦  K(K+l)q2^  ♦ 


A  A  A  r 

Z  - f ( Z  ,)[Kz*  . 
r+1  r+1  1  r+1 


C  R+ 1 ) q  j  1 
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V129 


c  lim  pr  E 


-(!)( 


383K 


P  -  (r-(2;*  -  r'(2:^!)]  - 


r(?i*r*l)] 


q.z  iK  In  z  +  1) 
n2  m  m 


*,.1«W 


fi'vi  *»  s.i  -  <K  in  h.i  *  nV 


=  lim  pr  E 


[r,,(2;zf)  -  r”(2;zJ5+  )]  2- 

P  .  _ JE _ - — •>  q  z  In  z  ♦ 

j  +  -  2  m2  m  m 

K  K 


W(V.iI  u„2 


-  r-K 
zr+l  Zr+1 


qxl> 


where  qx  =  r/n,  q2  *  (n-m)/n,  and  p  =  1  -  qj  -  q2  (m-O/n. 

As  n  +  «  with  q.  and  q  fixed,  zr+1  converges  in  probability 
1  z 


to  zr+1  Where  F(zr+1)  -  fQ 


f  (t) dt  =■  qx  and  z^  converges 


in  probability  to  z  where  1  -  F(z  )  e  f  t  f(t)dt  =  q2-  The 

m  m 

primes  indicate  differentiation  and  T(a;b)  is  the  incomplete 

Gamma  function  (Ref  7:559). 

Then  the  information  matrix,  in  terms  of  the  elements 

,  becomes 
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Let  a.,  be  the  elements  of  the  matrix  A,  which  is  the 
ij 

inverse  of  B,  then 
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(17) 


The  elements  of  the  variance-covariance  matrix  are  then 
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The  program  given  in  Appendix  C  finds  the  elements  of 
the  asymptotic  variance-covariance  matrix  by  solving  Eqs  (13) 
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(14),  and  (15)  to  find  the  elements  k^,  inverting  them  and 
converting  these  to  the  elements  of  the  variance-covarir  ce 
matrix  by  the  above  transformations. 

The  variance-covariance  matrix  is  solved  by  this  program 
for  non-integer  values  of  the  shape  and  scale  parameters 
since  with  the  location  parameter  known  the  equations  are 
regular  for  any  non-integer  value  of  the  parameter  estimates. 

A  flowchart  and  more  detailed  explanation  of  the  equations 
in  the  program  are  also  included  in  Appendix  C. 

Calculating  System  Confidence  Limits 

The  sample  values  of  system  reliability,  after  having 
been  generated,  are  ordered  and  thus  yield  the  simple  cumula¬ 
tive  distribution  of  the  system  reliability.  From  this 
distribution,  the  confidence  interval  and  limits  can  be 
approximated  at  any  level  of  confidence.  For  example,  in  an 
ordered  sample  of  99,  the  tenth  value  represents  the  lower 
limit  of  a  one  sided  confidence  interval  at  the  90%  confidence 
level  since  the  a  priori  probability  is  .90  that  a  new  random 
value  exceeds  the  tenth  order  statistic  of  a  sample  of  99. 
Thus,  it  is  assured  that  the  true  reliability  is  equal  to  or 
greater  than  this  value  with  a  confidence  of  .90. 

Verifying  the  Level  of  Confidence 

In  order  to  insure  that  the  confidence  level  is  accurate, 
a  program  was  written  to  generate  sample  times  to  failure 
from  the  Weibull  distribution.  Parameters  were  assumed  for 
each  of  the  component  distributions  3nd  various  size  sampler 
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of  tine  to  failure  were  generated  for  each  of  the  components. 
This  simulates  bench  testing  each  of  the  components  to  obtain 
a  sample  of  failure  times  for  each  component.  The  computer 
program  which  does  this  is  given,  with  an  explanation  of  the 
procedure,  in  Appendix  D. 

The  true  reliability  of  the  system  is  now  known,  sinc« 
the  component  distribution  parameters  are  known. 

The  failure  times  fox  each  component  were  fed  into  the 
Monte  Carlo  program  and  a  number  of  sample  values  of  the 
confidence  interval  at  the  required  level  were  generated. 

The  true  reliability  should  be  covered  by  a  percentage  of 
the  sample  confidence  intervals  equal  to  the  confidence 
level.  For  example,  if  10  one  sided  50%  confidence  intervals 
are  generated  and  the  true  system  reliability  is  .90,  a  graph 
of  the  intervals  might  appear  as  shown  in  Fig.  1  on  the  next 
page. 

In  this  case  half  the  intervals  cover  the  true  relia¬ 
bility  verifying  the  50%  confidence  level  (Ref  13:253). 

For  this  paper  confidence  intervals  were  generated  with 
each  interval  calculated  from  a  sample  (N)  of  99,  499,  999, 
and  2999  system  reliabilities.  The  results  are  given  in 
Chapter  III. 

Systems  Analyzed 

Two  systems  are  analyzed  in  this  paper.  The  first 
system  consists  of  8  components  connected  as  shown  in  Fig.  2 
for  which  the  required  mission  time  is  75  hours. 
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Fig.  2.  Schematic  of  System  A. 


The  components  have  the  following  characteristics: 
Component  >1. 

Failure  distribution:  Exponential 
Failure  test  sample  size:  100 
Sample  censoring:  top  6  ar.d  bottom  6  values 
Location  parameter:  0.0 
Component  #2. 

Failure  distribution:  Exponential 
Failure  test  sample  size:  75 
Sample  censoring:  Hone 


Location  parameter:  0.0 
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Component  #3. 

Failure  distribution:  Weibull 
Failure  test  sample  size:-  ISO 
Sample  censoring:  None 
Location  parameter:  0.0 
Component  #4. 

Failure  distribution:  Weibull 
Failure  test  sample  size:  200 
Sample  censoring:  None 
Location  parameter:  0.0 
Component  #5. 

Failure  distribution:  Exponential 
Failure  test  sample  size:  7S 
Sample  censoring:  None 
Location  parameter:  0.0 
Component  #6. 

Failure  distribution:  Weibull 
Failure  test  sample  size:  100 
Sample  censoring:  None 
Location  parameter:  0.0 
Component  #7. 

Failure  distribution:  Weibull 
Failure  test  sample  size:  200 
cample  censoring:  None 
Location  parameter:  0.0 
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Component  #8. 

Failure  distribution:  Weibull 
Failure  test  sample  size:  250 
Sample  censoring:  top  12  and  bottom  12  values 
Location  parameter:  0.0 

The  expression  for  the  reliability  of  this  '"stem  is 

Rs(t)  =  1  '  Vt)QbCt)  (20) 

where  Q.  -  1  -  R.(t) 

Ra(t)  =  [1  -  Q1(t)Q2Ct)][l  -  Q6(t)Q?(t)] 

Rb(t)  =  [1  -  Q3(t)Q4(t)Q5(t)]R8(t) 

The  second  system  (system  B)  has  the  configuration 


shown  in  Fig.  3  with  a  required  mission  time  of  100  hours. 


Fig.  3.  Schematic  of  System  B. 

The  components  of  system  B  have  the  following 
characteristics: 
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III.  Results 

Confidence  Limits 

It  is  known  that  the  error  associated  with  a  Monte  Carlo 
calculation  is  proportional  to  1//N  where  N  is  the  number  of 
simulations  (Ref  15:259).  In  this  case,  the  error  is  statis¬ 
tical,  that  is,  the  probable  error  is  proportional  to  1//N 
or  the  probability  is  high  that  the  approximate  solution 
does  not  deviate  from  the  true  solution  by  more  than  a  certain 
amount  (Ref  12:255).  The  use  of  a  high  speed  digital  computer 
should  decrease  the  amount  of  error  in  the  calculation  but 
the  machine  can  add  random  errors  of  its  own,  such  as  round¬ 
off  error  (Ref  9:11,12).  As  the  problem  becomes  more  complex, 
the  error  can  be  estimated  only  from  the  results  of  the 
computation . 

Tables  I  and  II,  on  the  following  page,  show  the  confi¬ 
dence  limits  determined  from  simulations  of  various  size  for 
systems  A  and  B,  As  the  simulation  size  increases  the 
accuracy  of  the  confidence  limits  at  a  given  confidence  level 
increases  and  more  closely  approximates  the  true  confidence 
limits.  Also  given  are  the  central  processing  times  for  the 
CDC  6600  computer  for  each  size  simulation.  No  trend  should 
be  noted  from  the  bounds  given  for  each  size  simulation  in 
these  tables  since  each  simulation  is  a  random  sample  of 
system  reliabilities  for  a  single  run.  The  only  proper  state¬ 
ment  that  can  be  made  is  that  the  confidence  limits  for  a 
simulation  of  size  2999  are  more  accurate  than  for  one  of 
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size  99.  Only  the  single  sided  confidence  limits  are  given 
in  these  tables  because  they  are  the  numbers  most  needed  in 
evaluating  the  reliability  of  a  new  system. 

There  is  also  an  error  associated  with  the  size  of  the 
sample  of  failure  times  for  each  component.  Kith  the  loca¬ 
tion  parameter  known,  the  sample  variances  and  covariances  of 
the  asymptotic  distribution  of  the  shape  and  scale  parameters 
are  different  from  the  variances  and  covariances  of  the  small- 
sample  non-asymptotic  distribution  of  the  parameters. 

Analyzing  this  error  is  difficult  but  Harter  and  Moore 
determined  the  error  to  be  relatively  small  if  the  location 
parameter  is  known,  even  for  samples  of  failure  times  as 
small  as  size  50  (Ref  7:562,  563). 

Lutton  compared  the  asymptotic  distribution  technique 
with  a  double  Monte  Carlo  simulation  developed  by  Moore  to 
approximate  the  system  reliability  confidence  interval, 

Moore’s  method  does  not  use  the  asymptotic  distribution  of 
the  parameters  but  uses  a  method  that  has  an  accuracy  that  is 
dependent  upon  the  number  of  simulations  made.  Lutton  found 
good  agreement  between  the  asymptotic  method  and  the  double 
Monte  Carlo  method  (Ref  11:21-26). 

Verification  of  Confidence  Level 

In  order  to  check  the  accuracy  of  the  confidence  level, 
the  proposed  method  was  repeated  several  times.  Here,  a 
single  repetition  of  the  method  is  called  a  run.  At  first, 
because  of  the  computer  time  needed  for  a  single  run  for 
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system  B,  100  runs  were  made  on  both  systems  A  and  B  for  a 
sample  of  99  system  reliabilities  (simulation  size  99).  In 
each  run  the  following  is  accomplished: 

1.  Sample  failure  times  for  each  component  are  gener- 

'*  ated  from  component  Weibull  distributions  with 

assumed  parameters. 

« 

2.  An  ordered  sample  of  99  system  reliabilities  is 
calculated  by  the  proposed  method. 

3.  A  counter  is  increased  by  one  if  the  true  reliability 
of  the  system  is  in  the  generated  confidence  interval 
at  any  given  confidence  level. 

After  all  the  runs  are  completed,  the  value  of  the 
counter  is  divided  by  the  total  number  of  runs.  This  gives 
the  percentage  of  runs  in  which  the  true  reliability  is  in 
the  confidence  interval  at  a  given  level;  as  in  the  example 
shown  in  Fig.  1.  The  results  of  this  procedure  are  given  in 
Table  III  on  the  next  page.  As  shown  in  this  table  the 
confidence  level  is  found  to  be  reasonably  accurate  for  a 
simulation  size  99,  which  is  the  least  accurate  interval 
because  of  the  error  associated  with  that  size  sample. 

Table  III  indicates  that  the  established  or  required 

* 

confidence  levels  in  the  proposed  Monte  Carlo  method  are  some¬ 
times  optimistic,  at  other  times  pessimistic.  It  must  be 
remembered  that  these  confidence  levels  were  checked  by  using 
only  100  runs  on  a  simulation  size  99  and  it  can  be  expected 
that  the  calculated  levels  would  approach  the  required  levels 
if  more  runs  were  possible.  Even  with  these  limitations, 
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Table  III  indicated  the  accuracy  of  the  proposed  method  is 
fairly  good  at  a  given  confidence  level. 

To  determine  the  effect  of  more  runs  on  the  calculated 
confidence  level.  System  B  confidence  levels  were  calculated 
for  300  and  700  runs.  This  was  done  only  for  System  B 
because  computer  time  for  each  run  is  not  excessive  for  this 
three-component  system.  The  results  are  given  in  Table  IV 
on  the  next  page.  In  most  cases,  if  the  calculated  level  was 
below  the  required  level,  the  calculated  level  approaches  the 
required  level  and  then  surpasses  it  when  more  runs  are  made. 
If  the  calculated  confidence  level  was  above  the  required 
level,  the  calculated  level  approaches  the  established  or 
required  level  when  more  runs  are  made.  In  almost  all  cases 
the  required  confidence  levels  for  this  system  tend  to  be 
conservative.  This  confirms  Lutton's  results,  ir  which  he 
found  the  confidence  intervals  for  the  asymptotic  distribution 
of  the  parameters  tended  to  be  wider  or  more  conservative  than 
in  Moore’s  double  Monte  Carlo  method.  In  other  words,  if 
the  results  in  Table  IV  apply  to  any  system  analyzed,  the 
confidence  limits  in  the  proposed  Monte  Carlo  method  tend 
(in  a  probabilistic  sense)  to  be  slightly  conservative  for 
both  one  and  two  sided  confidence  intervals.  This  is  a 
desirable  characteristic  from  a  practical  standpoint. 
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IV.  Conclusions  and  Recommendations 

Conclusions 

The  method  developed  is  extremely  flexible  and  gives 
good  accuracy  for  a  wide  number  of  system  configurations. 

The  only  inputs  needed  to  start  the  process  are  the  ordered 
censored  or  uncensored  samples  of  component  failure  times 
and,  along  with  the  known  value  of  the  component  location 
parameters,  any  kind  of  estimates  of  the  component  scale  and 
shape  parameters.  In  the  latter  case  a  guess  will  suffice 
since  these  estimates  are  used  only  co  start  the  method  of 
false  position..  An  increasing  vor  decreasing  failure  rate 
may  be  indicated  by  the  failure  times,  thus  supplying  a  guess 
for  the  shape  parameter.  Any  number  within  the  range  of 
failure  times  may  be  used  to  guess  the  scale  parameter. 

Bayesian  methods  or  Boolean  algebra  may  be  used  to 
express  system  reliability  in  terms  of  the  component  relia¬ 
bilities.  In  most  cases  this  can  be  done  with  a  single 
equation  which  is  inserted  in  the  main  program.  The  proposed 
method  can  be  applied  to  any  component  combination  and  uses 
relatively  little  computer  time  to  generate  a  large  sample 
from  which  confidence  limits  are  determined. 

The  accuracy  of  the  method  depends  mainly  on  the  number 
of  simulations  of  system  reliability.  It.  also  depends,  to  a 
lesser  extent,  on  the  size  of  the  sample  of  component  failure 
times.  As  the  size  of  either  of  these  increases  so  does  the 
accuracy  of  the  confidence  interval.  The  confidence  limits 


Recommendations 

It  is  recommended  that  further  work  be  done  on  the 
proposed  Monte  Carlo  method  in  the  following  areas: 

Types  of  Distributions .  The  method  should  be  extended 
to  include  other  distributions  such  as  the  Logistic  and  Gamma 
distributions.  Iterative  procedures  have  been  developed  to 
obtain  the  MLE  for  the  parameters  of  these  distributions  and 

could  be  incorporated  into  a  program  such  as  the  one  in  this 

\ 

paper . 

Error  Analysis .  It  is  also  recommended  that  a  wider 
variety  of  component  configurations  be  analyzed  to  try  to 
determine  if  there  is  a  direct  relationship  between  the 
number  and  configuration  of  the  components  and  the  error. 

In  other  words,  is  error  prediction  possible,  or  can  error 
be  determined  only  by  analysis  of  the  results? 

Analysis  of  systems  with  component  failure  time  samples 
of  less  than  ten  would  also  be  useful.  If  the  proposed 
method  remains  accurate  for  such  limited  component  testing, 
it  could  be  used  in  systems  which  have  components  with  a  long 
mean  time  to  failure  or  where  the  components  are  so  expensive 
that  only  a  limited  number  could  be  tested  to  destruction. 

Confidence  level  verification  as  given  in  Table  III  uses 
a  limited  number  of  computer  runs  on  simulation  size  of  99 
system  reliabilities.  It  is  recommended  that  further  examina 
tion  of  the  accuracy  of  the  confidence  level  for  simulations 
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of  499,  999,  and  2999  also  be  undertaken  to  check  the 
increase  in  accuracy  of  the  confidence  interval  for  these 
larger  size  simulations. 
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The  Parameter  Estimation  Subroutine 


T  |  <  —  1 T 

N=S AMFL7  c  I7-'  f:  N^G  - T )  , \  =  U  "  nn  LCSF  AT 

”1*"  Ic  'C  A  -  ft  f‘ :  7  r  ^  T-J-Tf,  T-  K-iOv’N 

5*1=1  IF  c  rt ..  -A  *  ^ '  T  ~  r  lrCC  tc  '?TIMAT-r 

*S2=.  T;  7  H F  .  A  *  V,rT" *r  <  73 

Sc;',~l  IF  ft-  \?r  ^  A  '•  V  '  T  “  ■'  ■<  T  ~  T7  ir  r  c.r  ;  »<  y n 

S  S  ?  =  ‘  IF  ir.:ATTru  :-t Tr  .  y  j--  •<jsr>vi, 

SF^-i  IF  L  C'A T 1 7>v  c-i;«v*T-v  £  y:  Tc  IF  -  *T  TMT^r 
TiTJsT-TM  T  A*  T  I C  ~r  $  ?••  =  (.£  (T=l 


r 

fTpP  .*  non 

C 

cm  =TMTT  I 

m  i  i  t  :  { 

V  \  L  Jc  ) 

>  ry  p 

r 

t.,htacm  =  :m 

t:ai.  ~v t 

c-  v 'L 

Vf 

r 

f  <  m.  )  =  i v:  i :  ?. 

L  _7TT-'r.r 

C  VTU1"* 

'’F  X 

r» 

yr  =  \  i ■>  .■  3  f 

-  uc.  7  t  c 

'o  '''a'"1':'  FfCM 

R~LCV» 

r 

v 

CUT"  LT 

r 

M  j  3  3 1  >  S  *  2 , 1  3 

'  71(17,  •"/(1)  »v'- 

r> 

--C 

»r  ,  r  p  r  _  t 

•  •'  3; 

r 

0  (J)=rSTIf'A~ 

'  1 F T r  •’  j-i 

•  *  1  T  V  N 

.  ‘  i  i.  '  . 

C 

(C  • 

v'f.r-.j  (  \ •  ’,  t  m 

i  '■‘T  n 

/  . 

C 

Tr;  F  T  A  (  J )  = r  '  1 

I ' '  -  *  li 

3-i 

r 

(C 

X'C-U,  «/!••• 

.)  '-r  '  >•  •'  7  7 

'y 

r<(  J3  =  FiTT> • 

T-  r  y  - 

•  : T  -  3  a  t  :  c  “ 

C 

r 


lr  4 
■j 


11  3 
115 


120 

123 

125 


(O'?  >'•*  l."  '  )  '-* 

( '*  Ay  *  jsjy  V  A  •_ !.:  ■  ]  r ~5.  -  ‘-ti  y  ry  «•  -►•ct  ONrO  Tc  r-') 

el  =K,AT’j‘ ii  lc-.  r.c  l : / - 1  * »  c  o  ,'th?ta  ( j) , r ►. "c  j> 

DT.xrv^ycN  T  (5  m  '  )  ,  T-1 17  M  •  ~  1 )  ,  -  v  n  ; )  ,  v  ( s6 )  .  Y  { ??  ) 

SFi-i, 

5*2=1, 

SS  3=  *  , 

IF  (4)  A  3  » 5  *  « <  V> 

5-M  =  M 

RE  A  3  3,  (T  (T) 
cO%v.\t  (  -  ?1 1  .  5 ) 

P- INI 

FPVA7M'I  ,lr'-  =  S7r”U“S«) 

P  r  1  \  T  CP 

FC  ■PAT  ( i  >t  ,  .«l  00.'  T  :f  u  ,  r  V  ,r  -J  '  7A  L‘  ,  5  v  ,r  Apr) 

0?lvT  l27,C(i)  ,TlltPMl)  ,•:•<{!) 

*0  PATCH 
P-T.t  123 
F 0  ? '•'  1 T  f  l;t 

EM  =  V 


r  £*/>,  . 

"  sQy. 


.•  / 


» i  ’  » 


u'i7'  -P.TiON*) 


31 

L  L  N  v  r  '  . 

E  H  R  =  v  r 

y  y  r.  r  v  ^  +  •< 

7  7 

V  -o 

NV:N*M*1 

00  ‘'4  7=N" 

FI  - T 

34 

LLN''  =  "L!.m-» 

IF  { r 

?  4 

no  /r  I  =  i  , 

HT  =  T 

F  L  N  w  =  •  ». \v-r LCG  ("T  ) 


i  3 


j 

U 

ii 

ll 


f  I 

y 


!  i 


i2J  fA3-.Os2A«w  AXarffr 


RAM/r,.A/72-2 


DO  ?w  J-l ,  t *?  “• 

IF  (J-l) 

JJ= J-l 
S'<  =  ‘  . 

SL  =  " « 

OC  r  I - v ^ ^ 

SK  =  SK4  (T(I)-C  (J  ])  >**F'<<J 
T.F  (?M) 

tht r ^  u)  -tk-t k  jj> 

GO  T0  0 

IF  <•'■>)  S':,  13,? 

TMET  1  ( J)  =  (  C-K  +  (<■  »•!-- 1  ) 

C** (1 ,/“<(uJ) ) 

GC  T0  <» 

x(1)=t^:ts  (  jj) 

hr**  Li,-.- 

L  L  = '.  - 1 
LP-L+i 
X  ( L  r )  =  X  ( L ) 

7;r<=  ( (T  hxF)-r  < j j>  )  n  ti  > 

Y  (L)  =-:>'(  Jj;  *  ( -  •--  ••  )  /v  ( 
c  +  *:  <<JJ><f(f  r.-F  v> *  (T  ( •  )-r  ( 

rj-'y-v 

IF  (Y  (L  )  )  :•*,  7  *•,='* 

Lc  =’  c - ' 

if  T l c  +  l )  r  5 ,  : ' 

LG  =  L'-4  ; 

Ir  (Lc-».)  :  ?,  J'  , 

X  (L°)  =  .r<v  (D 
GO  TO  ' 1 
X  (L'|rl,!*>'l  ) 

GO  TO  "'1 

I:  (Y  (L)*YC.L>  ) 

LL  =  LL -  * 

F-  C  T  0  r>  * 

X  {'."■)-<  (L  )  +v  (DM'  (i  )  -v  ( 
Ic  1 1-  ?  (X  v  L-)  -v  v  t  )  )  -!  .  c- 

roK tt»' 

T  «  -:  t  >  ( j )  -  y  ( i.  ~ ) 

EK(J) ==x ( JJ) 

IF  <SS’)  12,12,11 
n  0  O'7  I=“FF,“ 

SL*;UiLO'  <’(  T)-^(JJ)  ) 

X  ( 1  )  =  •  K  (  J ) 

Lc  =  - 

0  3  j  1  L  = 

F  L  *  -  .  • 

np  <;  Tr'^r,'1 

GlX-iLO  {  U  "G  (T  (  I  )-C  ( JJ> 
C**X  f l) 

LL  -I  -1 
LR  =  t  +1 
X  (Lr  ) =v  (L  ) 

7  R  K  =  (  (T  (*-•;  )  -•'  (  J  J)  )  /’l'rT 

Y  (L)  =  (0<'-f '•  <■)  *  i  1  . /V(L) 
f**y  (L)  +  ( ;:  S  - ' -^  )  *  (  -L0r  (’i:_ 
C-C  ( J  J  )  )  **  v  (Lj  .  r.*>  (  j ) * 
p  v  c  y  -7  (  _  t  r  \  /  ( 1  ,  -  '  y  ^  f  -  •»  V  i' ) 

IF  ( Y  ( L  )  ) 

L  5  =  L  c  -  1 

IF  ( L  r  +  L  )  ir^U-^L? 

L5  =  '.S*< 

IF  (L^-L) 

X  (L:i)  =  ,  G*y  fU) 

GC  TC  0 * 

V{l_F)r4.C*Y(L) 

go  io  ■;■ 

If  (  Y  (  L  )  *  y  ( lL  )  )  «-'■*<,  - ° r 

I.L-lL-1 

GO  7F  (*•» 


1-0  (  Jj)  )  **FF(  jj)  )  /ry) 


<(ji) 

i  )  4.  <(  JJ)  *C//Y<  L  )  *•>  (rK  (  JJ)  4-4  ,  ) 

)J))':-"-'(JJ)/X(u)'”^(rX(JJ)+l. 
"')  /  ( X  ( L  )  *  ( 1  #  -  :  X  F  ( -  7  3  K )  )  ) 


.1)  )/(>'  i'.J  *Y  ('»)  ) 
•  )  1 , 7 : , :  i 


)-'vl  OF.  (V-  =  TMJ)  )  )*(T  (I)-C(JJ)  ) 


-  I  J)  )  * *  *  { L  ) 

CO  T ( T  T 1 ( J) ) ) f-l-^LY/T^-T 1 ( J) 
T  J.  ( J )  4  -  ..LOG  (T  (  4)-C  (JJ)  ))  *  (T  (w) 
X  (  U  4 '  =  *  Z  F  rf*  (  j)i .  C 3  (  Z  SK  )  /  y  (  L )  ) 

) 


“s1'  ?:  >rf."V{f:»iiv5rif  VS^t'T/T?> 


GAM/MA/72-2 


X  (l_r)  -  X  (L)  4-y  (L)  ■•  (X  (t  )  -X  (LI  )  )  /  (Y  (U_>  -Y  (L>  ) 

IF  •  r  A""  ( y  ( L  =>1  -Y  l  L )  )  -i  .  i-  9)  *3  ?  >  c  u ,  =  1 
pontt^'j7 

FK  (  J)  =<  f'_r  ) 

C( j) =r. (  j j ) 

IF  /S'- )  /* 

IF  (l.-’S/o) )  :  ,75,7s 

IF  (  r  5  i  +  ~  c  ? 1  F  7  ,  f  7 ,  I  6 
x(i)--n(ji 

L  f=  ' 

DO  il  L  =  1  « c-  ' 

SK1=‘1. 

s=?=  . 

no  in  lrvrc,v 

SK1-SXT1  +  Cr  C)  -v  It  M  (r<l  J)  -  !•  > 

S-7  =  S^  +  1./  (T(I  )-*'  (L»  5 
LL  =  *  -1 
L  P  =  L  + 1 
X ( l c ) =v ft_) 

ZDX-((T(“'f)-v()  )  )  /•**<;  I  '  ( J)  )  **  "X  (  J  ) 

Y  ( L  )  =  (!«-  -  <(  J  4  :  v  < J  )  *  ( °  \  1  +  (  -  N  *  (  ’  < v  >  -  y  1  L )  ) 

("X  (  J )  -  ■  . )  )  /l ,,r  M  ( J)  **  "X  (  J )  -  I"  5* •"  <  ( J)  7^i<?ry  =  (_7  .  -.M 

0/  (  (T  (/'-")  -Y(D  )  *  ( 1  (-2  ?«') )  > 


37 

i c  ( v / l ) )  2'1, :*,/*■ 

L  S  =  1 3-1 

91 

IF  U_C*L>  N7  ,<*1,7'’ 

L3  =  l S+! 

91 

IF  {  L°-L  '  7  ,  ,  7  ' 

X  (LF1  =.  k*v  (i.) 

92 

r,C  TO  2'’ 

X  (  L  F  )  r  ,  -av  (•_)  +,-*T(; 

) 

73 

no  --  j 7 

IT  (V  (L)*Y  (LLn  7.  ,i 

- ,  7 1 

71 

L  L  =  L  L  -  1 

7? 

Q  p  'r  r\  -*  ' 

X  ( L  r  )  -  X  ( L  »  +  /  (  L  )  "  (ML. 

) -Y  ( 

2’ 

IF  !'»fc(x(L'.)-y(D)- 

i  •  '  ^ 

?o 

COi.T  T,.'JF 

?9 

C  (  J )  =  v  f  L  f  ) 

57 

GO  TC  7' 

C  (  J )  "  T  (  1 1 

25 

IF  (••■  )  i.j, 

35 

no  "3  1=1," 

IF  (r(J)H.’-t-T(Tl) 

67 

^.-.=w3  +  : 

6  3 

C  ( i )  =  T  (  1 ) 

65 

IF  («•')  •dF','9,7-' 

69 

s^=  • . 

Sl=  ' . 

F'O  ~6  I  =  “’"r,v 

S<rS<+  (’"(T)-''  (J>  ) 

<<  j) 

36 

S L  =  o L  +  •'■  L 1 r-  (I  (  T  )-F  /.!)  ) 

7ov  =  (  (T  C'^t  )-^{J)  )  /  ' 

»-  'T.. 

E.L  =  ‘  Lf  **  (  ■”  (  ■  • 

(•- 

C  (  r  <  (  j)  -  ;  , )  "l  -  (  h  ■»  ( -  )  MT  (*')  -0  ( J) )  X  ( J)  )  /  /  7r 

r  ( j ) w  „ ) )  *  .-.r-r-M.oc.  (i . -  .> -  (-■’ -x) ) 

F°  i v  T  ,0  U)  ,  P:T1  (  I)  ,  “<U>  ,rl 
?6  Ff’^STIlH  ,'.v,  ‘-i*./.,-!;  . ;  ) 


?.7  ir  1 a-- ;  <■:  («)-•:.  c  j  j  > ) -i -«•> 

?8  IF  TMJ  J)  )-t.  --'■*) 

jio  iF(d  ’••(-<(  j>  -  (JJ1  )  -1  ,  •  ,  1  2fc  ,  ^ 

30  CONTpM-: 

126  fnx-  T  iir 

1  *»  0  Cf)/'  M  (iIJ  ,  /  ,  If  ,  :  ■’.'-FI?'-/ 1  r  °  T  T  A  T  -  TSt> 

1  7 

U7  FT  -  T.T  {  1  I  ,/*<  ,  '  "l  TTF*  '.rat*  ,6X,  r,u -ru*  rr  ,r  Y. 

Cl.-iLCF  OF  IT,) 

pwi*  T  ,  F-  (J)  ,  T>  '  T  A  (  J '  5 :  <(  J )  PL 


37 


GAM/MA/72-2 


Appendix  B 

Generating  the  System  Reliabilities 

I. 

The  main  Fortran  program  SYSRE  generates  the  system 
reliabilities.  The  inputs  to  this  program  are  the  number 
of  components  in  the  system  (a  maximum  of  ten  as  presently 
dimensioned),  the  number  of  simulations  desired,  and  the 
mission  time.  Censoring  information  is  supplied  to  the 
program  along  with  the  known  value  of  the  location  parameter. 
The  failure  times  are  read  for  each  component  when  the  sub¬ 
routine  \PARES  is  called  and  the  true  value  of  the  location 
parameter  is  read  into  the  main  program  for  each  component. 
The  output  from  the  main  program  is  a  sample  of  ordered 
system  reliabilities  from  which  the  confidence  limits  can  be 
read  for  a  given  confidence  level.  The  flowchart  is  given 
in  Fig.  4  and  the  function  subprogram  which  is  used  to 
generate  pseudo-random  numbers  from  the  standard  normal 
distribution  is  explained  in  the  notes. 

Notes  (refer  to  Fig.  4): 

1.  The  ordered  component  failure  times  are  read  into 
subroutine  PARES  and  the  individual  component  MLE's 
of  the  shape  and  scale  parameters  are  calculated 
and  stored  each  time  the  subroutine  is  called. 

2.  When  the  covariances  are  equal  to  zero,  the  symmetr* 
matrix  to  vhe  1/2  power  is  simply  the  square  root  of 
the  non-zero  entries  of  the  2X2  matrix. 
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1/2 

3.  In  this  case  [V]  is  evaluated  by  using  the 
matrix  formula 

[V] 1/2  »  [S] ” 1 [A]n [S]  (22) 

where  [A]  and  [S]  are  matrices  that  diagonalize  [V], 
and  the  diagonal  elements  of  [A]  are  the  eigenvalues 
of  [V]  (Ref  3:244-251)  . 

4.  To  generate  samples  from  the  standard  normal  distri¬ 
bution  the  central  limit  theorem  was  used  where  the 
distribution 


>3) 


approaches  the  standard  normal  distribution  as  n 
If  t  *  x  -  1/2  where  x  is  a  pseudo-random  number  from 
the  uniform  distribution,  it  is  found  that  for  a 
sample  of  12  numbers 

12 

Y  -  l  x.  -  6  (24) 

i=  1  1 

has  a  distribution  that  approaches  the  standard 
normal  distribution  where  the  x^.  are  pseudo-random 
numbers  from  the  uniform  distribution  (Ref  9:145, 
146).  Function  subprogram  RANOR  accomplishes  the 
above  and  returns  the  numbers  to  SYSRE. 
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PROGRAM  SY5PP (INPUT, OUTPUT) 
r*¥***NCCMP  =  \'l;MPPR  CF  CCMrCNFNTS  IN’  SYSTEM 
C*****MJM=NUM”FP  OF  S I HUL fiTlC^S  DESIRED 
C*****TIFE  =  MSS3CN 

C*^^Mf»PEK  =  ESTIwS7E  CF  Sw  APE  PAOftMETcc  FCR  THAT  COMPONENT 
CJ,****ETH  =  ESTIf'A7E  0 p  p p A l. E  FARAMEIFR  FCR  THAT  COMPONENT 
C*»»v»xC  =  LCrATICN  FBc8f/ETE0  FOR  THAT  CCMp CKFNT 
C^**»»N«:AFP  =  SAFFLE  SIZE  0 F F C R E  CENSORING  ( N=r  3  C  C°  IPS*) 
C*****MASCV=CUANTTTY  N  HHE^r  f.-«  S  AM°L  E  ^  CENSORED  FRCH  APCVE 
C»<Mf**HRLCW  =  NUF?fR  OF  SAMPLES  '■'ENSCCED  FRCM  RrLCW 
PTMrNSTCN  S<2,2> ,V(1G,?»2)  ,R <5 £) » =c  (?C3 0 » * 

CNSAMP  (1C3)  .“A  )CV  (IOC)  ,WRLCK  (ICC)  ,  EEK{  13)  ,FTM10)  , 
CTXC(IC)  ,T)'(tO)  ,E<<1C)  ,TC <S5C ) ,E7 (55 C)  ,E«(55C) 
RS(1GC)=3. 

RP ( 5  0 1 ) =0  . 

RS(iuCQ)=C. 

rs  ( s  j  •:  n )  =  c . 

READ  E:,NCCvPtN'..'H  TIME 
r0  FORMAT  115  il1?  ie  2Z  »") 

PRINT  SP 

55  fCRMAT  (IF  ,1’WNC.  CP  CC^FC NENTS  , 5X  ,  »HNG .  OF  , 

C11WS  TNUI ATI0NS*CX*1?H FISSION  TIRE) 

PRINT  60,NCOM“,NU*1~IWF 
60  FCRWAT ( 1H  ,7X,I5,l?X,I5,l5X,F2n.5) 
on  13  2  X=i,NCOMF 
PRINT  63,  K 

6?  FCr,w,A'Sf  (IF  ,5° <*,*), //,1H  , IFFCC^FCN ENT  NlvrcR,ie) 

RE  AC  f  j,NP£MC(  ;<)  ,  vAFOV(K)  ,u-?LOW(K)  ,ETF(«<)  , 

CEEK(K) ,TVC  m 
70  FORMAT  (‘,Tc,3FQ.4) 

TC(1)=TXC  (*) 

ET  (1)=STP  (V) 

FKNtD-EE  KCO 

CALL  FARES  (NSAmc  (*)  ,KABOV  (<)  ,TC,E7, 
CFKN«.MrLCW(K),TurK),rN(K)) 

IF { Th ( < ) .EC.O.JPC  T C  RIP 
AA  =  vRLOK(1') 

Bpsp AnCV  (  X  ) 

CC=NSAVF  (K ) 

OlsAA/CC 

n?=(cc-en) /cc 

CALL  VACCV(Cl,C2,EK(<),CC»TH(K)»Vil,V12,V22) 

A  =  V11 

e=vi2 

C=V2? 

IF  (Vl2)!ir,l'10fll3 
100  V(K,1,1)=SCRT(V11) 

V(X,1,2)=0, 

V  <  K  ,  ?  ,  1 )  =  C  . 

V (<,2,2 )  =  S  CRT ( V2  2) 

V  ( 1 , 2 )  =  p . 

V  ( 2 , 1 )  =  p  * 

V ( 2 , 2 ) - SnFT (V22) 

GO  TC  112 

110  Y=. EFFORT (A*A-2.*A*C+C*C+4.*E*°) 

XI*. 5* (A+C)+Y 
X?=.5» (A+r )-y 

Rl  =  SCPf  (i.4(Xi-A>**2/8**2) 

R2=S0CT (1,4 (X2-0) **2/°**2) 

S(t,l)=l./Cl 
S(2,2)=l,/F2 
P(l,2) = (X2-R) / (n*R?) 

S(2,l)=(Xl-A)/(e*Rl) 

X1=SCRT (XI ) 

X2=SQRT (X?) 

D1=S  (1 , 1) *Xi 
D2=S(1,2>*>2 
03  =  p  <?,l)*vi 
C4=$(2,?)*X? 
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Appendix  C 


The  Variance-Covariance  Matrix 


Subroutine  VACOV  solves  Eqs  (13),  (14),  and  (15); 
inverts  the  result  and  multiplies  these  by  the  appropriate 
constants  as  shown  in  Eq  (17).  The  function  subprograms 
GAMI,  DGAMI ,  and  D2GAMI  compute  values  of  the  incomplete 
Gamma  function  and  its  derivatives  by  using  a  series  expan¬ 
sion  as  shown  in  Ref  6.  The  function  subprograms  GAM,  DGAM 
and  D2GAM  use  the  asymptotic  expansion  formula  to  calculate 
values  of  the  complete  Gamma  function  and  its  derivatives 
(Ref  4:257).  The  flowchart  for  subroutine  VACOV  is  given  in 
Fig.  5. 


Notes  (see  Fig.  5) : 

1.  vn  *  -  Kp  ♦  K(K+1)  [T  (2 ;  2^)  ]  +  K(K*l)q2£* 


P  .  tr-C2;S^J 

V22  ~  jT2  k2  +  q2"m  ln  Zm 


vu  ■  p  -  tr’(2;x*))  -  [r(2;^)]  - 


«Kri  „  0  -K, 

^2zm  *  *  k  *"  zm 


2.  vn  *  -  Kp  ♦  K  (K+ 1) 


v  -  p  .  r-(2.o 

V22  2  - 5 - 

K  K 


iflaSKSa* 
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SMI  ’OUT;-  7APCV  (  Ji,^ ?,'.•' ,  VU'I*  TW,?il  ,S12,$22> 

•>f  —  ”> 

J=  1 

TF  m.u'-'.l.JGC'  VP  1 
T F  (>.<:«.’  r  .  J  .  K-f  TC  1 
GO  "0  7 
1  fV  I v  T  ^ 

?.  Fn<vAl(l-i  ,  4  L  L  •!  <"•'  ‘>c  Ccrr‘) 

GO  TO  07 

7  rCMIf'J-  /X 

if  c::.:"...)  ^-o  -n  »  <v*X 

TL  Zc  =  r<*7F*,-  (.  <-j  ,  .)*  -<o(.7:«c<)  X^z/VV 

r-0  TO  q  X&x-X 

4  7^='.*  X 

y-\  ^ 

5  ir  XX. ' . > gc  tp  ^  NP 

7M=  v-^'.CG  ((*2)  )  **  ( 1  ,  .  */<'•*') 

GOTO  7 

7  7^*=4 .. : 

J= 1 

o2=i,r_oo 

6  Pal. 

if(u.i  i.nro  T~  it 

vn=-.  <*-*r<*  .  ) - g :  n :  )  +  '■«'*  (F'+i. . ) 

<**•  & ,,T  (7  *  <,  ?, '  ..>  /=■  <**?+,r,*7M**r<*3LC ; 

o  (  7  ■>  1  v?  7 

• '  -o-  ••*'  r»  :' 

r  (1  ,  •  .  -V*:  I_pr:  (  ?y)  ) 

r, r>  to  i  •> 

l  ^  v:t  =  -.  ';*r  *:-<*  r  •-<'*' .  ) 

V2?=T/rz^».  4.  ^  -  i  •:  (  -  ,  :  ,!/-^v7 

V  12  =  F-  ‘''''V('.  )  -  3  . 

17  IF  ('.  O,  i)  ;  TC  ^ 

vi i  = i <•( -  r  x*  r  '  3  r  )+,!*■-» f 

r  *rL7>?  rr^r^-v  <-  (  "  <  +  1 .  ‘  *  )  •  ~  1  )  'P  i*  *  ”> ) 

V27-Vi"-  <,  ■■■  '••  7  7- FL  7-^  A  ICG  (?->)  <  *'** 

V 1 2  =  V  t  2  +  0  o  *  ‘ T  ( 2  •  *  ’  ••  ,  w  . '  )  *  o  ’.  v  I  (  ?  ~  *  *  _  K  ,  2  .  ' )  -  r  1 

c  *:n  *=■•..  7  .*  :  **• i  r  ’)  -  <7* >  +  i. :  :> * 

POD  / 

?<+  0=VllvJ?2-VlwJf^'' 

Sll-VL'’/-1 

.si?=-vic‘/r 

S..1S0I  : 

S  2  °  =  V  _  1  /  0 

s i I = s r  imtf  -*.;i  /-‘.j- 

ci  >  =  » 

(  i  .  U”) 

P3I»‘T  2n 

26  Fc-r**’  (  i:-’  ,  ■’ttcriL';’  ]o..  7J  /jv-^cV  w -Uf  I  v  , // ,  1H  , 

C6X,  >v3',r>.7l-'':,,Ty1  ••HVr  7,/) 
p7['T  j-.,:.*, <i 
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p 7  ncNiir.-j- 
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>Vr^'CT'  ^-■■’7™;’.'!- 
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FU*wTI3rJ  f-  ft*!  (  Y ) 

7  =  v 

r«  • 

(z-°.'->  ?,:,r 

n  =  •;  -  A  I  0  r-  (  7  ) 

Z  =  ‘,«-i.  .  ■ 

o  o  y  p  J 

G:-  =  G  +(7-.‘  +  !*H0G(5«.‘»*3.i‘.lr- 

c  c’. !  .*Z  )-  '•  .  ‘/(7£-  .  -  7*  7***  )  *i  »t  ./ <12o..  . 

C-i.  '  .0  ■*"V¥7)+i  .  '  / 

is.7*;  v .  -  ' )  +  i.  -  “.  /  2) 

=:1’ 

F'J\~t:qn  c « 'i  y  (  y  >  >;^^Vsv 

•  v  ’  r  1 
,  j  •  .  :  1  •  j  .  / 

7  =  7  +  i  • 

r':'4rU(7.,=  ) /?*-:i'r  (7)  «i ..  -i. C ./ (I7.  .,*7**?) *i. r 

•'■/  ( i  ? .  ,  v  **  2**  h  )  -  4  i  .  /  ( Z'<: ■" .  .  ;*7**s)4-i.../(2i*  t  j ; *7**p> 
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Appendix  D 


Generatini 


>le  Component  Times 


Failure 


The  program  in  this  appendix  was  used  to  generate 
sample  failure  times  for  each  component.  If  the  points  in 
the  cumulative  distribution  function  are  represented  by 
random  numbers  in  the  uniform  distribution  over  the  interval 
0.0  to  1.0  (represented  by  Y)  then 


Y 


F(t)  =  1  -  exp 


(34) 


for  the  Weibull  distribution.  The  inverse  relation  is 


t  *  -0(-in(l-Y)]1/K  ♦  C 


(35) 


where  t  is  a  random  value  from  the  Weibull  distribution  and 
Y  is  a  pseudo-random  number  from  the  above  uniform  distribu¬ 
tion  (Ref  15:258).  Program  RSAMP  calculates  the  true  com¬ 
ponent  and  system  reliability  for  the  given  number  of  com¬ 
ponents  using  the  required  mission  time.  It  then  generates 
the  individual  component  failure  times  and  orders  them  so 
that  the  order  statistics  can  be  read  into  the  subroutine 
which  calculates  the  MLE  for  shape  and  scale  parameters. 

This  program  was  used  to  check  the  accuracy  of  the  computed 
confidence  interval  and  confidence  level  and  is  not  a  part  of 
the  method  proposed. 
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